Identification of a particle collision as a finite-time blowup in turbulence

We propose an Eulerian approach to investigate the motion of particles in turbulence under the assumption that the motion of particles remains smooth in space and time until a collision between particles occurs. When the first collision happens, particle velocity loses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document}C1 continuity, resulting in a finite-time blowup. The corresponding singularities in particle velocity gradient, particle number density, and particle vorticity for various Stokes numbers and gravity factors are numerically investigated for the first time in a simple two-dimensional Taylor-Green vortex flow, two-dimensional decaying turbulence, and three-dimensional isotropic turbulence. In addition to the critical Stokes number above which a collision begins to occur, the flow condition leading to collision is revealed; particles tend to collide in very thin shear layer constructed by two parallel same-signed vortical structures when Stokes number is above the critical one.


Identification of a particle collision as a finite-time blowup in turbulence Seulgi Lee 1 & Changhoon Lee 2*
We propose an Eulerian approach to investigate the motion of particles in turbulence under the assumption that the motion of particles remains smooth in space and time until a collision between particles occurs. When the first collision happens, particle velocity loses C 1 continuity, resulting in a finite-time blowup. The corresponding singularities in particle velocity gradient, particle number density, and particle vorticity for various Stokes numbers and gravity factors are numerically investigated for the first time in a simple two-dimensional Taylor-Green vortex flow, two-dimensional decaying turbulence, and three-dimensional isotropic turbulence. In addition to the critical Stokes number above which a collision begins to occur, the flow condition leading to collision is revealed; particles tend to collide in very thin shear layer constructed by two parallel same-signed vortical structures when Stokes number is above the critical one.
Particle-laden turbulence is frequently observed in daily life and industrial applications such as droplets settling in clouds and precipitation processes, particles sedimenting in river, aerosol pollutants in the air, fuel droplets in a diesel engine, and particles in a chemical reactor. Often, collisions between particles play a critical role in the determination of particle behavior. A good example is rain formation in which the collision rate between settling droplets is critically affected by background turbulence [1][2][3][4][5][6][7][8][9][10] .
For decades, the preferential concentration or clustering of inertial particles, which has been identified to be responsible for enhanced collision between particles, has been actively investigated using direct numerical simulation (DNS) [11][12][13][14][15][16][17][18][19] . The primary mechanism of the preferential concentration in the absence of gravity is the centrifuge effect of the inertial particles around rotating structures of turbulence. It is well characterized by the ratio of the particle response time scale τ p to the fluid time scale τ f , known as the Stokes number St = τ p /τ f . Especially when St ∼ 1 , preferential concentration is most pronounced in regions where fluid rotation is weak, but straining motion is not. The interaction between laden particles and turbulent coherent structures preferentially determines the accumulation pattern 20,21 . Moreover, when the settling motion of inertial particles due to gravity is dominant, a new type of preferential accumulation in columnar structures was observed 10,19,22,23 . Such a clustering naturally leads to more collision between particles. Although numerous studies on particle clustering revealed the physical mechanism and provided clues on how collision modeling can be improved, a full understanding of physical mechanism of collision lacks. Reference 24 attempted to investigate the collision mechanism tracking colliding finite-sized particles in direct numerical simulation of particle-laden isotropic turbulence. They showed that depending on Stokes number colliding pattern varies. However, an investigation of a collision event in the Lagrangian simulation of particles has limitation since it requires a luck to have a collision and thus it needs a large number of particles to be tracked to increase the chance of collision.
Equations of particle motion considered in this study are given by where r and v are the position and velocity of a particle, respectively. u is fluid velocity at the particle's position, which is a solution of Navier-Stokes equation and g is the gravitational acceleration vector. τ p is called the relaxation time of a small spherical particle defined by τ p = 2ρ p a 2 /9ρν with ρ p , ρ , ν and a denoting particle density, fluid density, fluid viscosity, and the particle radius. Equation (1) is valid under the point-particle approximation that the particle size is much smaller than the smallest flow length scale such as the Kolmogorov length scale when the density ratio of a particle relative to fluid is much larger than 1 so that all other forces such as the Basset www.nature.com/scientificreports/ history forces, the added mass and the fluid acceleration are negligible. It should be noted that when the density ratio is very large such as an order of 1000 (water droplet in the air), a small particle of order of 10 microns in cloud environment easily results in Stokes number of order one or larger. An inspiring characterization of particle collision was proposed by ref. 1 some time ago, who derived a formula for the collision rate of small inertial particles in turbulence using a form that traces random trajectories 25,26 . They introduced a particle velocity field v(x, t) from a solution for small Stokes number of equations of particle motion. Particle velocity field v(x, t) defining the velocity of a particle located at x at time t satisfies the equation ∂ t v + (v · ∇)v = (u − v)/τ p + g , which is equivalent to Eq. (1) under the assumption that the particle velocity field is a smooth function of space. Taking gradient of this equation yields, in a qualitative form, ∂ t σ + (v · ∇)σ + σ 2 = (s − σ )/τ p where σ = ∇v and s = ∇u . (This is made more precise in the next section.) When |σ | ≪ τ −1 p , it has a smooth evolution determined by σ (t) = t dt ′ exp[(t ′ − t)/τ p ]s(t ′ )/τ p . If |σ | > τ −1 p , however, σ 2 dominates and it may lead to an explosive evolution σ (t) ∝ (t 0 − t) −1 that produces singularity in v(x, t) in finite time. Although this estimate is based on the order-of-magnitude analysis, the singularity of the velocity gradient of a particle implies that the particle velocity develops blowup and becomes multivalued; in other word, a particle can have two different velocities at the same position, implying collision between two particles. Although compressibility of particle velocity or particle number density has been numerically investigated in turbulence using the Lagrangian method 6,[27][28][29] and the Eulerian model 30 , a clear identification of collision has not been carried out. Therefore, we propose a more rigorous Eulerian approach in which the particle velocity field v(x, t) is directly handled for an investigation of collision process of particles.

Results
Eulerian simulation of particle motion. We introduce a direct Eulerian simulation of particle motion within finite time. Equation (1) allows the following solution for the particle velocity, where w = gτ p is the settling velocity of a particle in still fluid. The particle velocity is determined by accumulated information of fluid velocity along a particle trajectory, which prevents further explicit expression. When τ p is small, i.e., a particle is small or light, however, the following approximation can be derived using the Taylor expansion of u(r(t ′ ), t ′ ) for small t ′ − t, Given that terms on the right-hand side of Eq. (3) are smooth functions of space, the particle velocity is uniquely determined by the position. Hence, the so-called flow of particle v(x, t) exists 10 . However, when τ p is not small, such an approximation is not possible because of the nonlocal nature of the solution, Eq. (2). As far as the particle velocity v(x, t) remains smooth, Eq. (1) can be written in the Eulerian form as Even when τ p is not small, Eq. (4) is valid under the condition that the particle velocity remains smooth because the left hand side of Eq. (4) corresponds simply to the Taylor expansion of the smooth particle velocity along a particle path. As will be shown later, the particle velocity remains smooth for a while when it was initially smooth until a blowup occurs. Then, the asymptotic behavior of the particle velocity leading to a blowup can be investigated using Eq. (4).
Taking the gradient of Eq. (4) yields an equation for σ ij (= ∂v i /∂x j ), where D/Dt(= ∂/∂t + v · ∇) denotes the material derivative along a particle path and s ij = ∂u i /∂x j . In the order-of-magnitude sense 10 , Eq. (5) reduces to where σ and s denote scalar variables representing σ ij and s ij , respectively. When τ p is small such that σ/τ p ≫ σ 2 , the quadratic term is small, yielding which remains finite since s is finite. When τ p is large such that σ/τ p ≪ σ 2 , however, the right-hand side of Eq. (6) is much smaller than the lefthand side, yielding a solution that might blow up in finite time, ∂v ∂t www.nature.com/scientificreports/ This kind of behavior leading to a blowup of the gradient of a solution is found in a solution to the inviscid Burgers equation when a shock develops 31 .
Reference 10 extended this analysis to further show that even when τ p is large, the finite-time blowup can be suppressed when gravity is strong or a nondimensional parameter called the Froude number defined by Fr ≡ η/(gτ 2 η ) is much smaller than one. Here, η and τ η are the Kolmogorov length and time scales. This implies that when gravity is sufficiently strong, the particle velocity may remain smooth. However, all these estimates are qualitative based on the order-of-magnitude analysis.
For a more rigorous investigation, Eq. (5), which is the tensorial form of Eq. (6), should be considered for analysis, indicating that the behavior of a solution cannot be simply estimated since the evolution of σ ij depends on the self-interaction between components of σ ij in a highly complex manner. It is also possible that quadratic nonlinearity does not always lead to a blowup, owing to the depletion of nonlinearity observed in the solution of the Euler equation [32][33][34] . Only a full numerical investigation of Eq. (4) can reveal the accurate behavior of a solution including the possibility of a blowup.

Three kinds of background flow considered.
To study a complete description of finite-time blowup of particle motion, we propose a direct simulation of Eulerian equation of particle velocity field, Eq. (4) for the given solution to the Navier-Stokes equation u(x, t) . In this study, we consider three types of the background flow in a cubic domain as described below. 2D decaying turbulence. Two-dimensional velocity field is given by u( The init i a l enst rophy sp e c t r um is sp e cif ie d by , where κ = |κ| with κ denoting the wavenumber vector. m and a m = (2m + 1) m+1 /2 m m! are the shape parameter and the normalization constant, respectively. We set u 0 = 10, m = 3, κ p = 10 and the detailed description is found in ref. 35 .
3D isotropic turbulence. In a periodic cubic domain [0, 2π] 3 , the three-dimensional velocity field u(x, t) is obtained from solving the continuity and Navier-Stokes equations, where p is the pressure, and f is a large-scale random force required to maintain stationary turbulence 15,36 . Parameters of f are chosen such that the Reynolds number based on the Taylor-scale, Re , is maintained at 14, which is intentionally set to be low to capture the blowup behavior well within the limited resolution. See "Methods" section for details of simulation methods and the definitions of the relevant nondimensional parameters.

Particle velocity fields near a collision.
To study the collision process of particles, we introduce a direct Eulerian simulation of particle motion within finite time (Eq. 4). Figure 1 shows the results of flow fields just before a blowup in the three background flows to provide the overall idea of a blowup visually. The region denoted by negatively large value of the divergence of particle velocity (Fig. 1a-c) or locally large value of the particle number density (Fig. 1d-f) corresponds to the place where a blowup is going to occur. As is wellknown 10,21,23 , particles with inertia tend to cluster in the region with weak fluid rotation but with strong straining motion. A blowup seems to occur in the region with strong compressive and weak stretching motion marked by a thin long region (Fig. 1), commonly for all three kinds of flow. Indeed, a blowup is easily identified in the Eulerian simulation of particle motion through an investigation of particle velocity field and particle number density distribution, although the simulation is valid only until the first blowup occurs since the simulation blows up numerically as well. In our investigation, we did not consider the size of a particle although the finite size shortens the collision time a little. www.nature.com/scientificreports/ Behavior of particle velocity divergence. A blowup can be identified in the behavior of the divergence of particle velocity (Fig. 1a-c). The evolution equation for the divergence of particle velocity is obtained through taking the trace of the evolution equation of particle velocity gradient (Eq. 5), since s ii = 0 due to incompressibility of fluid flow. The solution to Eq. (13) with the initial condition σ ii (0) = s ii = 0 can be expressed by indicating that the behavior of the quadratic nonlinear term along a particle path is playing an important role in the blowup process of σ ii . For small τ p , the integral can be approximated to yield where σ A ik is the skew-symmetric part of σ ik called the rotation rate tensor and σ S ik is the symmetric part of σ ik called the strain rate tensor. Given that σ A ik ≃ s A ik and σ S ik ≃ s S ik for small τ p , Eq. (15) suggests that an accumulation of particles quantified by negative divergence would occur in the place where the straining motion dominates the rotational motion, but this does not lead to a blowup.
When τ p is large, however, the exponential function decays very slowly in the integral of Eq. (14), allowing clearly indicating that a blowup could indeed occur through the quadratic nonlinearity of σ S ik since σ ii = tr(σ S ) in the particle accumulation process. Furthermore, σ ii = D j=1 j , where D is the dimension of space, and j is the eigenvalues of σ S ik with 1 > · · · > D . Therefore, the most negative eigenvalue D would blow up when a blowup of divergence occurs.
For an identification of a blowup, the temporal behavior of the L ∞ norm of | D | was numerically investigated for the three kinds of flow (Fig. 2a-f). When a blowup (see solid lines in Fig. 2a-f) occurs, D exhibits the following behavior commonly for all three kinds of flow, (4) obviously regularizes the solution by suppressing the divergence of the particle velocity since for small τ the particle velocity field tends to approach the fluid velocity which is divergence-free. When gravity is considered and thus particles settle, the blowup is delayed for the same St, and even a blowup which occurred in the absence of gravity is not observed (see cases with St = 0.4 for all three flows in Fig. 2a-f). This is consistent with the theoretical estimation by 10 in that as gravity becomes stronger, the flow of particles is observed in the wider range of St.
The critical time t c estimated from the behavior of 1/| D | as a function of St for various range of W (Fig. 2g-i).  Table 1. It is noteworthy that St c = 0.194 for W = 0 in 3D isotropic turbulence, while St c = 0.425 and 0.511 for W = 20 and 40, respectively, implying that the first collision between particles could happen at this critical Stokes number, although it might be rare.
Behavior of particle number density. When the divergence of particle velocity blows up, the particle number density also blows up (Fig. 1d-f). The conservation equation for particle number density n(x, t) reads which can be rewritten as www.nature.com/scientificreports/ Therefore, the solution for the uniformly distributed initial density n 0 can be expressed as , the solution of Eq. (20) can be obtained, or displaying the behavior similar to that of divergence or D . Numerical confirmation of this behavior will be given below.
Behavior of particle vorticity. If particle velocity is smooth, particle vorticity can be defined by ω p ≡ ∇ × v . From now on we discuss the evolution of ω p near a blowup. Figure 3 shows the fluid and particle vorticity fields just before a blowup for 2D turbulence and 3D isotropic turbulence. Near the position of a blowup marked by a black dot, fluid vorticity does not show any pronouncing behavior, while particle vorticity appears to be intensified for both 2D turbulence and 3D isotropic turbulence since ω p (x, 0) = ω(x, 0) with ω denoting fluid vorticity. Furthermore, a blowup seems to occur in a thin shear layer formed by two parallel vortices, commonly in 2D turbulence and 3D isotropic turbulence. In 2D Taylor-Green vortex flow, however, no such intensification of particle vorticity was observed since near the blowup, the flow field is irrotational due to the symmetric distribution (Fig. 1).
To understand this growing behavior of particle vorticity, we consider the evolution equation of ω p obtained by taking the curl of Eq. (4), The first and second terms on the right side imply the vortex stretching by straining motion and the vortex intensification by compression, which is usually absent in incompressible fields. Blowup happens in a thin vortex layer (Fig. 3), where converging motion is dominant over stretching motion. Therefore, the vortex stretching term can be neglected near the blowup. Then, multiplying ω p to both sides of Eq. (23) yields Near a blowup |∇ · v| ≫ τ −1 p and thus the last term on the right side of Eq. (24) is negligible. Then, for ∇ · v ∼ 1/(t − t c ) , the solution to Eq. (24) can be easily obtained, where C ω is a constant. This blowup behavior of ω p is similar to that of ∇ · v and n. Indeed, D , n and ω p show similar finite-time blowups with a proper constant (Eqs. 17,22,25). Numerically estimated constants for the three kinds of flow are listed for St = 0.6, 1 and 2 and various values of W in Table 2. All constants are found to be of order 1 for the range of St and W considered. Only C for 3D isotropic turbulence is around 2.5, which shows a little excursion from the values for 2D Taylor-Green flow and 2D decaying turbulence, while other constants Table 1. Critical Stokes number at various gravity factors for three kinds of flow. Critical Stokes number is obtained by fitting data (Fig. 2g-i) by t c /t f = a/ sinh(St − St c ) + b . The critical Stokes numbers for 3D isotropic turbulence are highlighted in bold face.
2D T-G vortex 2D decaying turbulence 3D isotropic turbulence  Fig. 4. This clearly confirms the near-blowup behavior of those variables for the three kinds of flow. Even with the relatively limited resolution for 3D isotropic turbulence, the blowup behavior is well captured although the behavior very near t c is untraceable further.    www.nature.com/scientificreports/ Condition for particle collision in 3D turbulence. Since a finite-time blowup implies a collision between particles, we can investigate the flow condition for which a collision between particles could occur. As we confirmed that a blowup occurs in a thin vortex layer of fluid and particle from Fig. 3, we examine more examples for universality of this scenario of collision. Since particle vorticity is usually unavailable, fluid vorticity distribution near a blowup for three different example 3D isotropic turbulence is provided in terms of the magnitude of vorticity in Fig. 5. A black dot indicates the location where it is on the verge of a blowup in each flow. It can be confirmed that a blowup or collision is likely to be occur in a pancake-type thin flat vortex structure induced by nearby vortical structures. It definitely serves as a qualitative condition for a collision. For a quantitative condition for a collision, we investigate the invariants of the fluid velocity-gradient tensor 37,38 . The fluid velocity gradient s ij (= ∂u i /∂x j ) can be decomposed into ω ω λ λ  www.nature.com/scientificreports/ with where s S ij and s A ij are the symmetric and the skew-symmetric parts of s ij , implying the strain-rate tensor and the rotation tensor, respectively. The eigenvalues S of s S ij , all real, satisfy the following characteristic equation, where P S , Q S and R S are the first, second and third tensor invariants, respectively, given by with where Eq. (30) is owing to the incompressibility of fluid and −4νQ S = ǫ , local dissipation rate. For s A ij , a similar characteristic equation holds with the corresponding invariants P A , Q A and R A , but P A = R A = 0 due to the skew-symmetricity of s A ij and which is positive definite. Therefore, A = 0, ±i|ω|. Distributions of joint probability density functions (PDF) between R S and Q S , and Q A and −Q S for 3D isotropic turbulence are demonstrated in Fig. 6 in which collision events observed in 10 different 3D isotropic turbulences are marked in colored dot. Two observations can be made: First, a collision tends to occur in the place of locally high strain but relatively low rotation, which is rarely found as shown in Fig. 6b. Second, most collision events are found closely to the right branch of the dashed line in the joint PDF between R S and Q S (Fig. 6a). The discriminant of the characteristic equation (Eq. 29), (27/4)R 2 S + Q 3 S , vanishes on both branches of the dashed line and sign( S 2 ) = sign(R S ) . Therefore, on the right branch of the dashed line, implying a straining field with one compressive direction and two equally stretching directions. All of these clearly suggest that a collision preferably occurs in a pancake-like vortex undergoing strong flattening.

Discussion
To investigate the near-collision behavior of particles in particle-laden turbulence, we studied an Eulerian form of equation of motion of particles (Eq. 4), which was derived under the assumption that the particle velocity is a smooth function of space. Using high-resolution simulations of the 2D Taylor-Green vortex flow, 2D decaying www.nature.com/scientificreports/ turbulence and 3D isotropic turbulence, we could accurately track a path leading to a blowup of the gradient of particle velocity, which was identified as a signature of particle collision. The temporal behavior of ∼ 1/(t c − t) predicted to occur in the gradient of particle velocity from an order of magnitude analysis by ref. 1 was found to hold in the divergence of particle velocity and particularly the most negative eigenvalue of the particle velocity gradient tensor when Stokes number of a particle is greater than the critical Stokes number which was also identified from our investigation. From a parametric study with Stokes number and gravity factor for the three different kinds of flow, we observed that as Stokes number increases, the blowup or collision occurs earlier and strong gravity delays the blowup.
When a blowup occurs, along with the divergence of particle velocity, particle number density also blows up in the similar manner, commonly for all three kinds of flow. Although our Eulerian simulation of particle motion is valid until the first blowup occurs somewhere, this direct calculation of particle number density could be very useful for multifractal modeling of uneven distribution of particles 10,20,21,39 .
Particle vorticity, which means spin velocity of a particle, exhibits similar blowup in both 2D decaying turbulence and 3D isotropic turbulence when the divergence of particle velocity blows up. From the equation for particle vorticity and simulation results, we identified that particle vorticity blows up through vortex intensification owing to compressive field of particle velocity. This is different from the finite-time blowup of vorticity of the solution to Euler equation through vortex stretching, which is possible only in 3D or axisymmetric flows [40][41][42][43][44][45][46][47][48] . It is interesting to notice from the equation for particle vorticity (Eq. 23) that the blowup due to compressive motion is suppressed in the incompressible Euler equation by pressure field, enforcing the divergence-free velocity field. In the particle motion, the Stokes drag force plays a similar role to pressure when Stokes number is below the critical Stokes number.
Finally, we investigated the flow condition most probable for particle collision in 3D isotropic turbulence. From various cases, we noticed that a thin flat vortical structure is the location where a particle collision is most likely to happen (Fig. 5). From a more detailed investigation of the invariants of the fluid velocity-gradient tensor, we confirmed that a collision preferably occurs in a pancake-like vortex under the compressive strain.
Our results of the near-collision behavior of particles are in principle valid for infinitesimally small particles. When the size of a particle is finite, however, the effect of size such as the lubrication forces between colliding two particles should be considered in the study of the near-collision behavior. Several collision models accounting for the lubrication forces before collision and the contact forces during collision have been proposed from particleresolved simulations based on the immersed boundary method [49][50][51][52][53] . However, given that the lubrication forces are non-negligible only when the gap between the colliding particles is much less than the radius of a particle, the near-collision behavior found in our study remains valid until the contact of two particles occurs as long as the particle size is much smaller than the flow length scale such as the Kolmogorov length scale.
We demonstrated from an Eulerian investigation of particle velocity that a blowup leading to particle collision can be accurately identified. This approach can be extended to an investigation of fractal structure of bubble motion in turbulence the behavior of which is different from that of heavy particles in that small bubbles hardly collide with each other due to small Stokes number.

Methods
Direct Eulerian simulations of particle-laden flows. Equations (4), (10), (11) and (12) are solved by a spectral method to maintain the spectral accuracy, and a time advancement is accomplished using a thirdorder Runge-Kutta method. Spatial resolution is kept at 4096 2 for the Taylor-Green vortex and 2D decaying turbulence and 256 3 for 3D isotropic turbulence. For the identification of the flow condition for a collision, an investigation of the invariants of the fluid velocity-gradient tensor was carried out using 128 3 resolution for 3D isotropic turbulence at Re = 70 . The Eulerian particle equation is solved in the frame moving at the settling velocity to avoid accumulation of numerical errors. In all cases, the initial condition for the particle velocity is v(x, 0) = u(x, 0) + w , which is smooth in space. The main nondimensional parameters are the Stokes number defined by St ≡ τ p /τ f and the normalized settling velocity W ≡ |w|/v f , where τ f and v f are the flow time and velocity scales, respectively. For the 2D Taylor-Green vortex and 2D decaying turbulence, τ f = 1/ω rms and v f = u rms where ω rms and u rms are root-mean-square vorticity and velocity of fluid, respectively. For the 3D isotropic turbulence, τ f = τ η and v f = η/τ η .

Condition for resolution.
To identify a blowup, we track large amplitude anomaly of the gradient of particle velocity, particularly the most negative eigenvalue of the symmetric part of the particle velocity gradient tensor since a blowup would occur when particles accumulate instead of creating vacuum. Figure 7 shows an example of a blowup in the 2D Taylor-Green vortex, in which the most negative eigenvalue 2 blows up in finite time. It clearly indicates 1/(t − t c ) behavior as t → t c (Fig. 7b) and as the resolution increases such behavior is better captured near t c (inset of Fig. 7a). Extrapolation of the linearly decreasing behavior of 1/| 2 | near t c also determines t c (Fig. 7a). Although an adoption of an adaptive mesh would yield better numerical evidence for a blowup 45,48 , the current resolution seems to be sufficient for the confirmation of the asymptotic scaling behavior as demonstrated in Fig. 4

Data availability
The data supporting the findings of this study are available from the corresponding author upon request.

Code availability
The simulation codes that have been used to produce the results of this study are available from the corresponding author upon request.